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ABSTRACT 


Based  on  theories  of  laminate  anisotropic  elasticity  and  interlaminar  frac¬ 
ture,  the  complete  solution  structure  associated  with  a  composite  delamination 
is  determined.  Fracture  mechanics  parameters  characterizing  the  interlaminar 
crack  behavior  are  defined  from  asymptotic  stress  solutions  for  delaminations 
with  different  crack-tip  deformation  configurations.  A  numerical  method  employ¬ 
ing  singular  finite  elements  is  developed  to  study  delaminations  in  fiber  com¬ 
posites  with  any  arbitrary  combinations  of  lamination,  material,  geometric,  and 
crack  variables.  The  special  finite  elements  include  the  exact  delamination  stress 
singularity  in  its  formulation.  The  method  is  shown  to  '>e  computationally  accurate 
and  efficient,  and  operationally  simple.  To  illustrate  the  basic  nature  of  com¬ 
posite  delamination,  solutions  are  shown  for  edge-delaminated  [6/-0/-8/6]  and 
[(:9) /(±8)  /90°/90<’  ]  graphite-epoxy  systems  under  uniform  axial  extension.  Three- 
dimensional  crack-tip  stress  intensity  factors,  associated  energy  release  rates, 
and  delamination  crack-closure  are  determined  for  each  individual  case.  The  basic 
mechanics  and  mechanisms  of  composite  delamination  are  studied,  and  fundamental 
characteristics  unique  to  recently  proposed  tests  for  interlaminar  fracture  tough¬ 
ness  of  fiber  composite  laminates  are  examined.  Influences  of  lamination,  geome¬ 
tric,  and  crack  variables  on  the  delamination  behavior  are  investigated. 
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In  an  associate  paper  (11,  the  fundamental  nature  of  stress  singularity 
and  asymptotic  solution  fields  associated  with  a  delamination  between 
dissimilar  anisotropic  fiber  composites  have  been  studied.  Lekhnitskii 's 
complex-variable  stress  potenc.'^ls  [2]  in  conjunction  with  an  eigenfunction 
expansion  method  have  been  used  in  the  formulation  and  establishment  of  the 
general  solution.  The  eigenvalues,  especially  the  ones  which  characterize  the 
strength  of  stress  singularity,  for  delaminations  with  different  local 
traction  boundary  conditions  near  the  crack  tip  have  been  determined.  Of 
particular  Interest  are  the  asymptotic  deformation  and  stress  governing  the 
composite  delamination  fracture.  Since  logarithmic  terms  are  absent  in  the 
solutions  for  the  delamination  stress  and  deformation  [1],  the  general 
structure  of  the  complete  solutions  consists  of  only  a  power-type 
eigenfunction  series  of  both  singular  and  higher-order  terms.  The  set  of 
unknown  constants  in  the  eigenfunction  series  solution  for  a  delamination 
problem  can  only  be  determined  by  solving  the  complete  boundary-value  problem 
with  a  full  consideration  of  overall  composite  geometry,  lamination  and 
material  variables,  remote  boundary  conditions,  and  end  loading  conditions. 

In  a  finite-dimensional  fiber  composite  laminate  with  simple  lamination 
variables  and  crack  geometry,  for  example,  a  symmetric  angle-ply  [8/-0/-0/0J 
composite  containing  edge  delaminations  along  the  8  and  -6  ply  interface ,  the 
complete  laminate  elasticity  solution  can  be  determined  in  an  explicit  form  by 
the  use  of  different  analytical  methods  such  as  the  boundary-collocation 
method  (3).  However,  for  a  composite  laminate  having  more  than  four  plies, 
the  aforementioned  collocation  method  is  not  applicable,  and  a  more  general 
and  versatile  method  of  approach  to  the  problem  is  needed.  The  situation 
could  become  extremely  complex  for  delaminations  in  a  composite  containing  a 
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large  number  of  plies  wlch  different  fiber  orientations  and  laminar 
thicknesses.  The  almost  unlimited  number  of  variables  In  a  general  composite 
delamlnation  problem  requires  the  development  of  an  advanced  analytical  method 
which  can  not  only  deal  with  the  local  singular  behavior  of  the  delamlnation 
but  also  take  Into  account  various  combinations  of  lamination,  material, 
geometric,  and  crack  variables.  Numerical  methods  such  as  finite  element 
methods  are  considered  to  be  most  attractive  because  of  their  versatility  In 
handling  mechanics  problems  with  complex  structural  geometry  and  material 
properties.  In  this  paper,  the  currently  developed  laminate  elasticity 
solution  Is  Incorporated  in  the  formulation  of  a  set  of  conforming  finite 
elements  with  singular  derivative  fields  [4].  The  special  crack-tip  elements 
are  shown  later  to  be  particularly  suitable  for  modeling  the  composite 
delamlnation  problem.  The  use  of  laminate  elasticity  solutions  derived  from 
the  associated  paper  [1|  permits  the  Inclusion  of  exact  delamlnation  stress 
singularities  In  finite  element  formulation,  leading  to  extremely  accurate  and 
efficient  numerical  solutions  for  studying  the  fundamental  behavior  of 
composite  delamlnation  with  complex  lamination  variables  and  geometric 
parameters. 

Specific  objectives  of  this  paper  are  to:  (1)  establish  complete 
solution  structures  for  different  delamlnation  configurations  to  serve  as  a 
basis  of  formulating  an  advanced  numerical  method ,  (2)  construct  special  crack 
elements  of  various  kinds  to  model  the  composite  delamlnation,  (3)  properly 
define  Interlaminar  fracture  mechanics  parameters,  e.g.,  stress  Intensity 
factors  and  energy  release  rates,  for  general  composite  laminates,  (4)  examine 
the  fundamental  behavior  and  associated  characteristics  of  the  composite 
delamlnation,  and  (5)  assess  influences  of  lamination,  geometric  and  crack 
variables  on  the  delamlnation  response  and  composite  failure  modes. 
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In  chc  n«xc  section,  the  general  structure  of  complete  solutions  for 
delaainatlon  stress  end  displaceaent  fields  Is  given.  Fracture  mechanics 
parameters  In  terms  of  Interlaminar  crack-tip  stress  Intensity  factors  and 
strain  energy  release  rates  are  defined  for  various  deformation  modes  of 
composite  delamlnatlon.  Special  six-node,  quasi  three-dimensional  crack 
elements  for  the  present  composite  delamlnatlon  problem  are  Introduced  In 
Secclon  3.  Formulations  of  the  singular  elements  and  adjacent  nonsingular 
eight-node  Isoparametric  elements  are  briefly  outlined.  Solution  strategy  and 
computational  scheme  for  the  delaainatlon  problem,  especially  in  the  case  that 
crack-surface  closure  (l.e.  the  contact  problem)  occurs,  are  discussed. 
Computational  methods  for  evaluating  stress  intensity  factors  and  energy 
release  rates  by  using  the  singular  finite  elements  are  given  in  detail.  In 
Section  4,  solution  accuracy  and  convergence  are  studied  to  demonstrate  the 
efficiency  and  effectiveness  of  the  present  approach.  The  fundamental 
behavior  and  unique  characteristics  of  composite  delamlnatlon  are  examined  In 
Section  5.  Two  graphite-epoxy  laminate  systems,  l.e.,  symmetric  angle-ply 
[3/-8/-8/8J  and  symmetric  l(±8)/(±9)/90#/90*l8  composites,  with  delamlnatlons 
emanating  from  laminate  edges  are  studied.  Delamlnatlon  crack-tip  deformation 
and  fracture  mechanics  parameters  are  determined  for  each  case.  Influences  of 
fiber  orientation,  ply  thickness,  and  crack  size  on  delamlnatlon  failure 
mechanics  and  mechanisms  are  investigated  also. 
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2*1  Complete  Solutions  for  Stress  sod  Displacement 

As  discussed  In  the  associated  paper  (1],  Che  eigenvalues  4„  can  be 
determined  by  solving  the  transcendental  characteristic  equation  Eq.  20  in 
(1|.  Depending  upon  the  local  crack-tip  boundary  conditions  and  the  interface 
continuity  conditions,  the  6a  have  the  following  values: 

(1)  delaai nation  with  open  crack  surfaces 


fn  -V2 


4  •  (  n 


(n  -  V2  )  *  it 


(single  root), 
(triple  roots), 
(single  root). 


(n  •  0,1,2,...*); 


(11)  delanlnatlon  with  closed  crack  surfaces  In  frictionless  contact 


fn  -V2 

n 


(double  roots), 
(quadruple  roots). 


(n  •  0,1,2,...*); 


(ill)  delanlnatlon  with  closed  crack  surfaces  in  frictional  contact 


S 

n 


-l/2 


-  0 


(single  root), 
(single  root), 
(quadruple  roots). 


(n  *  0,1,2,...*). 


Once  the  values  of  4n  are  determined,  the  relationship  among  C^'s  can  be  found 
and  the  complete  solutions  for  displacement  and  stress  can  be  established  In 
explicit  forms.  For  example,  for  a  delanlnatlon  with  homogeneous  local 
boundary  conditions  in  a  composite  laminate  subjected  to  planar  loading,  the 
complete  solutions  for  stress  and  displacement  components  have  the  following 
expressions : 
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orrs.""J_  ?; [\ 

OF  POOR  QUALITY 


i  •  L(ci»  J. <iu  ik<“  vr 1/2 1  *  vr 1/2 " 


n-0  k-1 


\^c2nbto  *3«bin  »  "^Ik#  +  *c4nbi„  I  ^Ik2^ 


♦  <„!>  *lt2J'V2*lTl  ♦  Vj,  C  #,kC1/2+1''l»  +  ».l 

(i  -  1,2, 4, 5, 6),  (1) 


(k*3>n!'(JM)kritn  !•/<*»  Vi) 


♦  J,  < « *«•  «r< j«>k*r > '  ♦  =«„bta)l-lr<j+«kzfl'<"*»'> 


*  b™  !wzrV2Mr'<-l'2«»>i 


k-1 


+  c 


‘"'■'J,  6b”>  r(j*3)k2rl/2Mr/<»*i/2«',||,i  *  %j  «> 


where  b£*^  ere  known  eigenvectors  corresponding  to  the  unknowns  can  for  each 

given  6n,  and  oot  and  uQj  are  known  quantities  froa  particular  solutions  for 

(a) 

each  Individual  case.  We  note  that  the  b^n  are  found  to  be  the  saae  for 

all  Vs  because  of  the  lnvolveaent  of  the  ten  ein<*(4  -  *  and  -*)  In  the 
n 

D(&n)  aatrlx.  The  constants  A1Jt  and  are  defined  as 
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Alk"  uk»  A2k“  l*  A4k"  “V  A5k"  wkV  A6k“  "V 

rik"  1 »  r2k“  V  P3k“  uk*  r4k“  pk»  r5k“  ‘•k*  r6k“  ck*  (3) 

utatrt  Ujj,  pjj,  qjj,  tjj  are  related  to  ply  stiffness  matrix  and  can  be  found 
In  (2,3).  The  expressions  for  stresses  and  dlaplaeeaents  for  delaal nation 
with  closed  crack  surfacee,  l.e.,  nonvanishing  local  tractloa  boundary 
conditions,  can  be  easily  determined  to  have  forms  similar  to  Eqs.  1  and  2  but 
with  slight  modifications  owing  to  different  numbers  of  algebraic  multiplicity 
of  the  eigenvalues  Involved. 

2.2  Del  ami  nation  Stress  Intensity  Factors  and  iMrgy  la  lease  lates 

The  stress  and  displacement  fields  for  a  delaninatlon  are  shown  la  the 
preceding  section  to  possess  general  form  of  Eqs.  1  and  2  with  unknown 
constants  to  be  determined.  A  proper  analytical  or  numerical  method  with 
the  aid  of  global  lanlnate  boundary  conditions  and  renote  loading  conditions 
Is  required  to  determine  the  detailed  solution  for  the  complete  boundary-value 
problem.  Since  the  Interlaminar  fracture  is  controlled  by  local  stress  and 
deformation,  the  asymptotic  solution  Is  of  primary  importance  and  Interest  In 
understanding  the  near-fleld  behavior  and  fracture  phenomenon.  The  asymptotic 
solution  Is  recognised  to  be  singular  In  nature  and  governed  by  the 
delamination  stress  singularities,  which  have  been  obtained  In  detailed  In 
(1).  As  pointed  out  In  the  associated  paper  (1),  the  singular  eigenvalues 
depend  upon  the  local  delaninatlon  configuration;  thus,  distinct  structures  of 
asymptotic  solutions  are  obtained  for  different  crack-tip  deformations. 

In  general,  the  asymptotic  solution  for  s  delaninatlon  stress  field  can 
be  written  In  the  following  form: 


l 


V 


7 


t 
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J-l 


If 


k-1 


<J>  t  J 
lk  *k 


+ 

ll(k+3) 


* 


(1  ■  li2|3|i(>|i}i 


(4) 


where  and  have  their  origin  located  at  tba  da laminae Ion  tip. 

Iha  fjj[*  and  f^+3j  are  function#  of  lanlna  aatarlal  conatanta,  lamination 
variables,  and  loading  eondltlona.  tha  lntagar  t  la  tha  total  nuabar  of 
eigenvalues  6j  which  satisfy  tha  conatralnt  condition, 

-I  <  «a(«j J  <  0  (5) 

For  tha  convenience  of  furthar  developments,  Iq.  4  may  ba  rawrlttan  aa 


•  l  •i1)<r,d;  4  ), 
1  j-l  1  J 


(6) 


where  a^  la  tha  jth  coaponant  (corraapoodlng  to  tha  eigenvalue  6j)  of  tha 
aaynptotlc  atraaa  ot  in  tha  polar  coordinate  ayatan. 

In  view  of  tha  asymptotic  aolutlon  structuraa  given  in  Eq.  4  and  thoae 
given  In  [1,3,6],  It  la  poaelble  to  defloe.  In  tha  context  of  Interlaminar 
fracture  aachaolca,  tha  dalaml nation  atraaa  lntenalty  factors  and  strain 
energy  release  rates  In  a  aanner  consistent  with  those  for  a  honogeneoua  crack 
and  for  the  refined  modal  of  an  Interface  crack  between  dissimilar  Isotropic 
madia  [7,8).  For  example.  In  tha  case  of  a  closed  da lamination  with  crack 
aurfacaa  In  frlctlonless  contact,  tba  stress  Intensity  factors  can  be 
Introduced  by  considering  the  erack-tlp  Interlaminar  stresses  oj,  and  o4 
(i.e.,  Oy,  Tgy,  tyX,  or  o^,  rr^,  T(^>  along  the  ply  Interface  ♦  -  0  as 
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•  lin_  ^-2*r  a,( r,0) 

(7a) 

r*0  2 

•  llm*  /2«r  o.(r ,0) 

(7b) 

r*0  * 

•  llm.  /2tr  o4(r,0) 

(7c) 

r*0 

where  t  •  0“  and  0*  denote  th«  positions  laflnlteenally  behind  and  ahead  of 
cha  crack  dp,  respectively,  because  the  normal  traction  11  finite  ahead  of 
a  doaed  crack  tip  where  shear  tractions  og  and  04  are  singular. 

The  strain  energy  release  rate,  C,  and  Its  components  Ct  (1  •  1,2,3)  may 
be  evaluated  by  using  Irvin's  virtual  crack  extension  concept  (9)  as 
C  •  Cj  ♦  Cj  ♦  Cj 


-  lie  ^  Mo^r.OJlu^tfa-r,.)  -  u<**l)<fa-r,-s)| 
fa  *0  0 

♦  o^r.OMu^tfa-r.*)  -  u(I*fl><fa-r,-v>| 

♦  o4(r,0)(u£B,(fa-t,»)  -  tt^Bt',)(dB-r,-s)|)dr,  (8) 


where  fa  la  the  length  of  virtual  crack  extension.  The  interlaminar  stresses, 
«2,  Og,  and  04,  in  89.  8  cas  be  obtained  from  the  asymptotic  stress  field  such 
as  89.  4.  The  corresponding  displacements  are  also  those  of  the  asymptotic 
field  e9uatlone  discussed  la  the  previous  section,  la  terms  of  the 
<!e lamination  stress  Intensity  solutions  Kj,  the  C  and  Cj  for  a  closed  crack 
can  be  shown  to  have  a  simple  expression  as 

c  "  ®2  ♦  ®3  “  7  <A2*Il  *  A3Km>  <*> 


In  89.  9,  C|  is  Identically  aero  because  of  the  dlsplacemsat  continuity  across 
the  closed  delaml nation  surfaces. 


tt(n)^rtt)  _  B(m*-I)(r>_tj  ,  0 


-c  <  r  <  0, 


(10s) 


\ 

9 

and  *2  **4  AJ  lr*  svaluaced  fro*  IS]  with 

“i*5  "  “2**1*  “  (10b) 

UjB)  -  u<~l)  -  k^l\io~r\  (10c) 

»»  t  *  a*. 

to*  a  do lamination  with  extreaely  snail  crack  closure,  tha  slnpllfled 
model  by  taking  tha  Halting  can*  of  a  partially  closed  crack  discussed  In 
(5,61  la  osad.  Tha  across  Intensity  solutions  Kj [  and  Kjjj  are  the  saae  as 
thoaa  glean  In  Eqs.  7(b)  and  7(c),  but  tha  Kg  la  daflnad  it  r  ♦  0+  as 

IL  -  lia.  fui  o,(r,0).  (11) 

1  r*0  1 

dad  tha  corresponding  strain  «n*rgy  ralaaaa  rates  C  and  Cj,  chan,  have  the 
fora. 

c  •  c,  ♦  C2  ♦  C3  -  i  <*,*2  *  *2*1!  ♦  A,«4n).  C 12) 

to  the  ease  that  a  fully  opened  da lamination  la  assuaad  and  the 
eigenfunction  expansion  series  Is  used  for  Che  ssyaptotlc  solution,  the 
de leal  nation  stress  Intensity  factors  are  Introduced  as  (3| 

3  ... 

%  -  11a  I  /5t  r  3  Sj  (r,0;4  ),  (13a) 

r*0  J«1  J 

3  - i ... 

-  li*+  I  /T7  r  3  e£J,(r,0;«.>,  (13b) 

T*0  J-l  3 

3  ... 

<4n  *  11%  I  SR  r  J  s;J,(t.0;«  ). 
r*0  J»1  1 


( 13c) 
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We  note  chat  In  chis  case  exacc  Integration  in  Eq.  8  can  be  carried  out  and 
the  strain  energy  release  rate,  G,  can  be  determined  easily.  However, 
Individual  components  G^  can  not  be  separated  explicitly  In  the  integration 
because  of  the  complex  mathematical  structure  of  the  asymptotic  solutions. 
Other  methods  such  as  the  hybrid  singular  finite  element  analysis  [10,11]  are 
needed  to  determine  the  values  of  individual  . 
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3.  COMPUTATION  OF  FIELD  SOLUTIONS  AND  FRACTURE  MECHANICS  PARAMETERS 


We  have  thus  far  obtained  the  complete  solutions  for  stress  and 
displacement  and  the  fracture  mechanics  parameters  from  the  asymptotic 
solutions  for  each  individual  case.  The  unknown  constants  c^  as  well  as 
and  are  to  be  determined  by  solving  the  complete  boundary-value  problem, 
involving  remote  loading  conditions  and  global  geometric  and  lamination 
parameters.  Now  consider  a  finite-width  composite  laminate  containing 
delaminations  under  mechanical  loading  as  shown  in  Fig.  1  in  [1].  As 
mentioned  in  Section  1,  several  numerical  methods  could  be  used  to  determine 
the  complete  solution  for  the  delamination  problem  in  a  finite  dimensional 
laminate.  Owing  to  the  singular  nature  of  the  interlaminar  crack,  the  complex 
structural  geometry,  and  numerous  lamination  parameters  involved,  the  special 
numerical  method  employing  recently  introduced  conforming  singular  finite 
elements  is  an  attractive  approach  to  the  current  problem.  In  this  approach, 
Che  exact  delamlnatlon  stress  singularities  can  be  included  in  the  formulation 
of  the  special  elements.  Thus,  the  unknown  constants  associated  with  the 
asymptotic  solution  and  the  fracture  mechanics  parameters  governing  the 
delamlnatlon  behavior  can  be  evaluated  conveniently  with  a  high  degree  of 
accuracy  and  a  fast  rate  of  solution  convergence. 

3.1  Singular  Crack-Tip  Elements  and  Surrounding  Nonsingular  Elements 

In  this  study,  we  generalize  the  formulation  and  concepts  of  the 
conforming  singular  elements  originally  introduced  for  homogeneous  isotropic 
elasticity  problems  [4]  to  the  present  quasi  three-dimensional,  anisotropic 
composite  delamlnatlon  problem.  Formulation  of  the  singular  elements  is  based 
on  selection  of  shape  functions  and  their  derivatives  containing  the  exact 
eigenvalues  which  meet  the  constraint  condition  Eq.  5  of  stress  singularity 
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derived  from  the  aforementioned  eigenfunction  analysis  for  the  composite 
delaaination.  Detailed  discussion  of  the  nature  of  the  singular  eleaents  has 
been  given  in  (A].  Only  relationships  relevant  to  the  current  developaent  are 
given  here. 

Consider  a  six-node  triangular  eleaent  witc  three  degrees  of  freedoa  per 
node  and  stress  singularity  at  node  1  as  shown  In  Fig.  1(a).  The  element  has 
general  rectilinear  anisotropic  material  properties  with  elastic  compliance 
Sjj.  By  a  proper  transformation,  any  point  in  the  element  defined  in  global 
Cartesian  coordinates  (x,y,z)  can  be  referred  to  both  local  polar  coordinates 
(r,t,z)  and  triangular  coordinates  (p,£,z)  with  the  origin  located  at  node 
1.  Within  the  eleaent,  the  displacement  components  u  are  related  to  the  nodal 
displacements  q  by  the  Interpolation  (or  shape)  function  as 


u  *  U(x,y)  +  U  with  U  **  N  q,  (14a, b) 

—  «  .o  «  „s  w 

where  uq  are  known  quantities  resulted  from  applied  loading,  and 


T 

u  ■ 

{  ul» 

V  u3}-  “o  "  {uol*  Uo2*  Uo3}« 

(15a) 

(1) 

_(  !) 

J2)  „(2)  J2)  „(6)  „(6)  -(6)1 

ql  ■*  q2  *  q3  »•••••*  »  ^2  *  q3 

(15b) 

2  ’ 

q3  ' 

Ns  -  (P,  CJ  «j>. 

(15c) 

in  which  the  superscript T  denotes  transpose  of  the  associated  column  vector,  and 
the  nuaber  in  the  superscript  parentheses  in  Eq.  15(b)  refers  to  the 
associated  nodal  nuaber  of  q.  Explicit  expressions  of  the  shape 
function  N#,  Involving  proper  eigenvalues  <5j  determined  for  each  Individual 
delaaination  problem  from  [1]  and  local  coordinates,  are  given  in  Appendix 
1.  Equation  14(b)  can  be  written  in  a  more  explicit  form  as  follows: 

V1 

U  -  {p  J  MU)  +  L(p.O)  q. 


(16) 
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where  L(p,€)  is  Che  nonsingular  polynomial  pare  of  H (p,E;6,),  Independent  of 
-  J 

Che  number  and  location  of  nodes  along  Che  2-3  edge  of  Che  element.  Using  Che 
minimum  potential  energy  theorem  and  following  the  same  procedure  as  the 
conventional  displacement-based  finite  element  formulation,  we  can  construct 
the  element  stiffness  matrix  kg  for  the  crack-tip  element  as 

kg  -  JJ  £  C  Bs  dA,  (17) 

A 

s 

where  C  Is  the  material  stiffness  matrix,  and  B  has  the  form 

.  .s 

6+1 

Sg  -  3(p  J  M  +  L).  (18) 

In  Eq.  *8,  the  3  is  a  matrix  differential  operator.  We  remark  that  the  shape 
function  of  the  singular  element  is  chosen  such  that  the  element  conforms  with 
a  nonsingular  quadratic  element  matched  through  the  common  element  boundary 
(l.e.,  along  edge  2-3),  and  with  singular  elements  of  the  same  formulation 
through  boundaries  1-2  and  1-3. 

The  surrounding  nonsingular  elements  used  in  this  study  are  quasi  three- 
dimensional,  eight-node  isoparametric  elements  [Fig.  1(b)]  with  24  degrees  of 
freedom  (three  D.O.F.'s  per  node).  Formulation  of  the  element  stiffness 
matrix  for  the  nonsingular  element  has  been  given  in  detail  in  [12,13].  The 
element  stiffness  matrix  kf  for  the  adjacent  element  can  be  shown  as 

k  -  //  8^  C  B  dA,  (19) 

.r  A  -r  .  .r 

where 

Bf  -  3  Nf.  (20) 

The  shape  function  Nr  for  an  eight -node,  quasi  three-dimensional  isoparametric 
element  has  standard  quadratic  expressions  which  can  be  found  in  [12,13]. 
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3.2  Solution  Procedure  end  Iteration  Scheme  for  Crack  Cloanre 

In  the  finite  elenent  discretization  of  the  delaminated  composite,  the 
singular  elements  are  placed  In  a  ring  fora  with  the  delamlnatlon  tip  as  the 
common  node  (Fig.  2).  The  element  stiffness  matrices  k  and  k  are  true 
stiffness  matrices  relating  unknown  nodal  displacements  to  nodal  force 
vectors.  The  standard  procedure  of  the  matrix-displacement  method  [14]  can  be 
used  to  assemble  the  global  stiffness  matrix  K  and  loading  vector  Q  leading  to 
the  relationship. 


K  q  -  Q  +  Q0. 


(21) 


where  Qq  denotes  the  additional  nodal  force  resulted  from  the  applied  strain 
e0,  and  the  assemblage  may  be  expressed  symbolically  by 


K 

r* 


“s 

l 

i-1 


(22a) 


"S 

Q-  l 

~  i-1 


(22b) 


in  which  n8  and  nr  are  the  total  numbers  of  singular  and  nonsingular  elements, 
respectively. 

In  actual  n>  aerlcal  computation  of  a  delamination  problem,  a  finite 
length  of  crack  closure  is  assumed  first.  The  problem  now  becomes  an  elastic 
contact  problem  because  a  part  of  the  local  boundary  conditions  is  not  known 
and  needs  to  be  determined  from  the  solution.  Specifically,  the  following 
continuity  conditions  are  required  along  the  closed  portion  of  the 
delamlnatlon,  -c  <  r  <  0: 
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(a2<r,w)]  -  0,  (23a) 

[u2(r ,*)]  -  0,  (23b) 

where  che  brackec  [  ]  denotes  the  jump  of  the  associated  quantity  across  the 
closed  crack  surfaces,  e.g.,  uj^^r  ,-*)  -  u£m^(r,»)  *  0.  Thus,  we  solve 
Eq.  21  in  terms  of  unknown  pressures  along  the  closed  crack  surface  and  use 
Eq.  23(a)  to  determine  the  contact  stress  on  this  surface.  To  enforce  the 
conditions  Eqs.  23(a)  and  23(b),  the  solution  technique  proposed  by 
Francavilla  and  Zienkiewicz  [15]  for  an  elastic  contact  problem  is  employed. 
The  numerical  procedure  Involves  an  iteration  scheme  to  determine  crack 
closure  length  and  contact  stress  along  the  delamination  surface.  If  the 
solution  is  admissible,  the  contact  stress  so  obtained  must  be  in  compression 
and  the  displacement  field  should  have  no  overlapping  or  interpenetration 
outside  the  contact  region.  In  the  case  that  crack  closure  length  is  found  to 
be  extremely  small,  say,  less  than  the  order  of  10  5  '•  10  6  times  crack 
length,  the  delamination  is  then  assumed  to  be  open,  and  the  simplified  model 
discussed  in  [1,5]  and  in  the  preceding  sections  in  this  paper  is  used.  For 
the  case  of  a  delamination  with  finite-length  crack  closure,  the  detailed 
Iterative  algorithm  for  evaluating  crack  closure  and  contact  stress  is  given 
in  Appendix  2. 

3.3  Computation  of  De  lamination  Stress  Intensity  Factors  and  Energy  Be  lease 

Bates 

A a  mentioned  in  the  preceding  sections,  the  stress  intensity  factors  and 
energy  release  rates  for  a  delamination  are  evaluated  from  the  asymptotic 
solution  of  interlaminar  stresses  <Jg  and  (or  oy,  x^  and  xyz)  and  the 
displacements  ut  along  the  plane  of  the  crack.  For  a  finite  dimensional 
composite  laminate  containing  delamlnatlons ,  the  asymptotic  stress  and 
displacement  can  be  conveniently  determined  by  che  aforementioned  singular 
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finite  element  method.  Along  the  delemination  plane  5  -  CQi  the  near-field 
stress  and  displacement  fields  are  approximated  by  using  Eq.  16  and  its 
derivatives  as 

flap6  R(5o;5>  q,  (24a) 

u  =  p6+1  M(£  )  q  +  u  .  (24b) 

.  _  o  _  _o 

where  R(€  ;6)  is  a  matrix  of  derivatives  of  the  shape  function  M  and  ply 
>o  _ 

elastic  constants,  and  the  p  and  50  are  related  to  the  global  coordinates  by  a 
simple  transformation  given  in  Appendix  1. 

For  a  partially  closed  delamination  or  a  delamination  with  a  very  small 
size  of  crack-tip  closure  for  which  the  simplified  model  with  inverse  square- 
root  singularity  is  used,  the  stress  intensity  factors  and  strain  energy 
release  rates  can  be  evaluated  easily  from  the  singular  finite  element 
results.  Taking  a  and  u  along  the  delemination  crack  plane 
Co  -  0  (i.e.,  -  0)  ,  we  can  write  the  asymptotic  Interlaminar  stresses  and 

displacements  along  the  interface  in  simple  expressions  as 

o  *  A.  r  ^  (i-2,  k-1;  i-6,  k-2,  and 

1  1-4,  k-3),  (25a) 

u.  *  B.  r  +  u  .  (j-1,  k-2;  J-2,  k-1,  and 

J  j-3,  k-3),  (25b) 

where  Afc  and  are  obtained  from  the  corresponding  components  of  R(€q;6)  and  M(£q)  in 
Eqs.  24(a)  and  24(b)  by  setting  5Q  -  0,  and  5  -  -  Vjj  •  Thus,  the  delamination 
stress  intensity  factors  can  be  easily  determined  by 

K1-/IFa1  (i«I,II,III),  (26) 

The  energy  release  rates  can  be  determined  in  a  manner  similar  to  that  for 
through  Eq.  8  as 
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G  -  Gi  +  C2  +  g3 


11a  1 

5a +0  25a 


•  £o  Jr  I  I.W")'**-  «’•> 

0  k«l 

The  Cera  uQj  In  Eq.  25(b)  is  not  Included,  because  It  does  noC  result  in  any 
contributions  to  G^  and  G  after  Integration.  Integration  of  the  singular 
Integral  Eq.  27(a)  can  be  carried  out  explicitly  without  difficulty.  Then  the 
strain  energy  release  rates  have  the  fora 


/**{  l  (\  r'^2)[Bk(«a-r)  ^JJdr 


0  k-1 


6a  3 


r  5a -ry^2 


G  “  J  Gi  "  T<A1BI  +  V2  +  SV* 


We  reaark  here  that  each  term  in  Eq.  27(b)  corresponds  to  the  individual  G^ 
components  and  that  for  a  delaaination  with  finite  crack  closure,  the  first 
term  in  Eq.  27(b)  is  Identically  zero,  l.e.,  Gj  -  0,  because  of  the  continuity 
of  dlsplaceaent  across  the  crack  surface. 
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4.  SOLUTION  ACCUIACT  AID  COMVERGEICE 

To  demonstrate  Che  accuracy  and  efficiency  of  Che  presenc  method  of 
approach,  a  symmetric  angle-ply  {45®/-45®/-45®/45® J  graphite-epoxy  laminaCe 
containing  edge  delaainadons  along  the  45*/-45*  ply  lnterace  Is  considered. 
For  simplicity  and  without  loss  of  generality,  the  composite  Is  assumed  to  be 
subjected  to  uniform  axial  strain  ez  »  e0  along  the  z-axis  and  have  a  geometry 
of  b/h  -  8  and  hj  -  t»2  "  h  ■  0.25  Inch  and  delaminations  of  length  a  ■  0.25 
Inch  emanating  from  the  edges.  Elastic  ply  properties  of  unidirectional  high- 
modulus  graphite-epoxy  identical  to  those  of  Eq.  28  In  [1]  are  used.  Owing  to 
the  geometric  and  lamination  symmetry,  only  a  quarter  of  the  cross  sectional 
area  needs  to  be  considered.  In  finite  element  discretization  of  the 
continuum,  twelve  special  crack-tip  elements  of  identical  size  and  shape  (Fig. 
2)  are  used  to  model  the  near-field  response  of  the  composite  delamination. 

The  crack-tip  elements  are  embedded  in  the  mesh  of  eight-node,  quasi  three- 
dimensional  isoparametric  elements.  Local  and  overall  mesh  arrangements  for 
the  finite  element  analysis  of  the  composite  delamlnatlon  are  shown  in  Fig. 

3.  To  study  the  accuracy  and  convergence  of  field  solutions,  the  mesh  near 
the  crack  tip  is  continuously  refined  by  halving  the  lengths  of  equal  sides  of 
the  singular  elements  (e.g.,  the  sides  0E  and  OF  in  Fig.  3). 

Using  the  computational  method  and  the  solution  scheme  discussed  in  the 
preceding  sections,  numerical  results  are  obtained  for  the  convergence 
study.  Significantly  global  crack  closure  with  c/a  -  0.34  is  found  for  the 
delamination  in  the  (45®/-45®/-45®/45®)  graphite-epoxy  under  the  uniform  axial 
extension,  resulting  in  a  negative  Kj  and  identically  vanishing  G^.  {The 
detailed  nature  of  delamination  closure  and  related  problems  will  be  discussed 
in  the  next  section.)  In  Tables  1  and  2,  delamination  stress  intensity 
factors  and  energy  release  rates  associated  with  the  finite-element  mesh 
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refinement  are  shorn  for  systematic  reduction  in  crack-tip  element  size 
(Ah/h)  and  Increase  in  the  total  number  of  elements.  It  can  be  seen  clearly 
that  stable  and  converged  solutions  for  all  and  Gt  are  obtainable,  as  the 
crack-tip  element  size  becomes  smaller  than  0. Ih  and  the  associated  total 
number  of  elements  exceeds  100.  As  anticipated,  the  delamination  stress 
intensity  factors  (Table  1)  are  sensitive  to  the  size  of  the  crack-tip 
elements  owing  to  the  localized  nature  of  the  singular  domain.  However,  the 
strain  energy  release  rates  are  relatively  Insensitive  to  the  element  size  and 
mesh  refinement  (Table  2)  because  the  G^  are  related  to  the  global  structural 
response  of  the  delaminated  composite.  Furthermore,  the  crack-closure  length 
is  also  found  to  be  insensitive  to  the  mesh  refinement  in  this  case. 

To  assess  the  accuracy  of  the  current  solutions  for  the  composite 
delamination  with  global  crack  closure  is  not  trivial  because  no  analytical 
and  numerical  solutions  are  available  in  the  literature  for  comparison.  The 
only  analytical  study,  which  may  be  used  as  a  reference,  deals  with  the 
delamination  problem  by  assuming  the  crack  surface  being  fully  open  and 
employing  an  eigenfunction  expansion  method  with  the  aid  of  a  boundary 
collocation  technique  (3|.  In  Fig.  4,  stress  intensity  solutions  determined 
by  the  present  singular  finite  element  method  including  the  crack-closure 
consideration  are  presented  as  a  function  of  delamination  length  in  the  (45*/ 
-45*/-45*/45*J  graphite-epoxy.  The  dominant  stress  Intensity  factor  Kjjj 
determined  for  the  partially  closed  composite  de1.  ;  nation  by  the  present 
approach  is  about  6  -  7Z  higher  than  the  value  for  an  open  crack  from  the 
boundary  collocation  results.  However,  solutions  for  Kj  obtained  by  using  the 
two  distinct  models  differ  from  each  other  both  in  sign  and  in  magnitude.  In 
the  current  singular  finite  element  analysis  of  the  delamination  problem,  Kj 
is  found  always  to  be  negative  owing  to  the  aforementioned  crack  closure. 
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whereas  Kj  determined  for  an  aaauaed  open  crack  is  positive.  Values  of  Itjj 

and  G2  determined  by  the  two  approached  are  both  negligibly  smell.  The  strain 

energy  release  rates  obtained  for  the  delaalnatlon  are  given  in  Pig.  S.  For 

the  partially  closed  delaalnatlon  determined  by  the  present  approach,  Gj  Is 

6  2 

Identically  zero,  whereas  the  value  of  G^IO  cq)  from  the  assumed  open  crack 
model  Is  of  the  order  of  10~^ — a  value  much  smaller  than  that  of  G3,  the 
dominant  component  of  energy  release  rates.  Despite  the  different  nature  at 
the  crack  tip  In  these  two  approaches,  the  values  of  G3  differ  only  slightly 
from  each  other. 

We  remark  that  though  mathematically  rigorous,  the  solution  obtained 
using  the  eigenfunction  expansion  method  with  the  assumption  of  open  crack 
surfaces  Is  physically  Inadmissible  because  severe  Interpenetration  of 
dissimilar  materials  Is  found  in  the  oscillatory  displacement  solution  near 
the  crack  tip.  Thus,  the  present  singular  finite  element  approach  with  a 
partially  closed  crack  consideration  provides  physically  meaningful  solutions 
with  a  high  degree  of  efficiency  and  accuracy  for  general  composite 
delamlnatlon  problems. 
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5.  THB  VODAMOrrAL  BEHATIOK  UD  OUIACTWSTICS  Of  COMPOS IT*  PHJMHJMATIOM 

la  Che  preceding  eeccloae,  we  have  introduced  the  coeplete  eolution 
structure  end  defined  the  governing  fracture  aechanics  parameters  for  the 
delaaination  problea.  Special  finite  eleaenta  with  singular  derivatives  have 
been  foraulated  to  aodel  the  near-field  response  of  a  delaaination  in  a  finite 
diaensional  coaposlte  laainate.  Solution  convergence  and  accuracy  have  been 
afflraed  in  teras  of  coeputatlonal  peraaeters  (e.g.,  the  else  of  crack-tip 
eleaenta,  the  degree  of  finite  eleaent  discretisation,  etc.).  Having 
established  these  basic  aechanics  theories  and  the  nuaerlcal  aethod,  we  now 
proceed  to  study  the  fundaaantal  behavior  and  characteristics  of  delaainatloo 
in  fiber  coaposites  with  genaral  lamination  and  geoaatrlc  variables. 

Two  fiber  coaposlte  aaterlal  ays tea*  are  examined  in  this  section:  (l) 
symmetric  angle-ply  [±0]g  graphite-epoxy  laminates  with  edge  delaalnatlons 
between  0  and  -6  plies,  and  (2)  symmetric  l(±9)/(i$)/90*/90r}B  graphite-epoxy 
with  edge  delaalnatlons  between  -0  and  90*  plies.  The  symmetric  angle-ply 
coaposlte  laainate  system  is  selected  because  several  unique  delaaination 
characteristics  are  observed,  which  can  be  used  to  illustrate  most  clearly  the 
basic  Interlaminar  fracture  aechanics  and  failure  nodes.  More  laportsntly, 
soas  of  the  aost  fundamental  nature  of  delaaination  fracture  determined  froa 
the  presently  Introduced  physically  admissible  aodel  and  mathematically 
rigorous  solutions  are  not  observable  In  the  previously  obtained  solutions 
which  contain  the  Inadmissible  oscillatory  stress  and  deformation  [3,5].  The 
l(*0)/(±«)/9O */90*]g  graphite-epoxy  system  Is  studied  also  because  this 
lamination  system,  especially  the  one  with  l(*30,)/(*30*)/90*/95irIs  fiber 
orientations,  is  currently  belag  considered  for  use  In  the  evaluation  of 
interlaminar  fracture  toughness  of  composite  materials  under  static  and  cyclic 
loading  [16,17],  The  ply  elastic  properties  of  hlgh-aodulus  unidirectional 
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graphl te-epoxy  given  la  Eq.  28  of  8ef.  (1]  art  used  la  all  computations.  la 
Che  (8/-8/-8/8)  composites,  cac  following  laalaacloo  and  gvonttrlc  parameters 
art  employed:  hi  *  h2  "  h  “  0.25  Inch  and  b/h  -  8,  whtrtaa  la  the 
((*8)/(e6)/90*/90*)g  composite  laminates,  actual  dlatoalona  of  unldlrtccloaal 
graphite-epoxy  used  la  laboratory  experiments  (16,17)  art  taken,  l.e*, 
hi  -  hi  ■  ...  »  hi i  •  0.007*  Inch  and  b  -  0.75  inch.  Moreover,  la  all  eaaae 
atudltd  la  thla  paper  the  composite  laalaatea  art  coaaldertd  to  be  aubjtcted 
to  uniform  axial  exttnaloa  c{  ■  cQ  along  the  t  direction. 

Since  the  lacerlaainar  crack -tip  deformation  aad  fracture  meehanlsma  are 
governed  by  the  aayaptotlc  field  aolutlona,  ue  shell  examine  the  fundamental 
behavior  and  aaaoclaced  characterlatlca  of  compoalte  delaml nation  la  terma  of 
Interlaminar  fracture  mechanlca  parametera,  l.e.,  crack-tip  atreaa 
latenaltlea,  atraln  energy  releaae  ratea,  and  crack-earface  cloture. 

5.1  Influence  of  Piker  Orientation 

The  behavior  of  a  delaml  net Ion  la  the  relatively  almple  angle-ply  |*8)g 
laminate*  la  algnlflcantly  Influenced  by  the  fiber  orientation  8.  Reaulta 
obtained  by  ualag  the  Iterative  solution  achent  In  Section  3.2  reveal  that  the 
delaml  net Ion  always  poaaeaae*  a  finite-length  crack-tip  eloaure.  Assuming 
that  the  crack  surface  la  In  frlctlonleaa  contact,  we  find  that  the  closure 
length  la  global  In  general.  For  example,  In  the  (±8),  graphite-epoxy  with 
edge  delaml  nations  of  length  a  -  0.25  Inch,  significant  crack  eloaure  la  found 
In  each  case  studied  (Table  3).  The  cootact-sooe  site,  c/a,  varies  from 
approximately  two-tenths  to  more  than  one-third  of  total  delaml nation  length 
for  8  ranging  between  IS*  to  60*.  The  crack-tip  closure  reaulta  In  a  negative 
opening-mode  stress  Intensity  factor  <  0  and  an  Identically  vanishing 
energy  releaae  rate  Cj  •  0,  as  shown  In  Tables  6  and  5.  Thus,  the 
de lamination  behavior  In  the  (8/-8/-8/8)  composites  is  apparently  governed  by 
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Interlaminar  shear  scrcsse*  end  deformations.  In  Table*  4  and  5,  Che  out-of- 
plane  tearing-mode  (node  111)  screas  intensity  factor  and  associated 

energy  release  rate  Gj  are  clearly  the  dominant  ones  controlling  the 
delamlnaclon  fracture.  (The  values  of  Kjn  and  Cj  are  several  orders  of 
magnitude  higher  than  those  of  and  Gj.)  Also,  the  composite  with  0  *  lb' 
possesses  the  highest  Kjjj  and  C}  among  all  8's  studied.  These  results 
suggest  that  la  symmetric  angle-ply  composites,  delamination  Initiation  and 
growth  are  more  intimately  associated  with  local  Interlaminar  shear  then 
transverse  normal  stress. 

Ue  reiterate  chat  the  solutions  obtained  hers  by  the  use  of  the  present 
InterlasUnar  crack -closure  model  and  the  rigorous  mathematical  (combined  theo¬ 
retical  and  numerical)  approach  are  physically  ad^sslble  and  meaningful,  as 
contrary  to  the  previously  obtained  Inadmissible,  oscillatory  solutions  from 
an  eigenfunction  expansion  approach  with  Che  assumption  of  a  fully  open  crack. 

In  the  more  complex  graphite-epoxy  laminates  with  [(*6)/(*6)/90*/90*)s 
ply  orientations,  delamlnaclons  are  always  observed  to  occur  between  the  -0 
and  90*  plies  |16,I7).  Starting  with  the  partially  closed  crack  model  and 
using  the  Iterative  solution  scheme,  we  find  that  the  closure  length  of  a 
de lamination  Is  less  than  10~*  inch.  Thus,  for  this  extremely  small  crack 
closure  the  slapllfled  model  with  an  Inverse  square-rooc  stress  singularity 
discussed  In  Section  1.3  of  (1|  and  la  (S,6|  Is  employed.  In  Figs,  b  and  7, 
variations  of  stress  intensities  and  energy  release  rates  with  fiber 
orientation  9  are  shown  for  a  del ami net  ion  of  length  a  •  0.5b.  The  crack  tip 
Is  apparently  governed  by  the  opening-  and  inplane  shear lng-a»de  stresses;  the 
and  are  one  order  of  magnitude  higher  than  for  all  8*s  studied. 

The  value  of  Gj  is  vanishingly  small  (n  general,  end  C|  and  G2  are  three 
orders  of  magnitude  higher  then  Gj  because  of  the  extremely  small  value  of 
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A3.  Note  that  values  of  Gt ,  G2  and  G3  depend  not  only  on  hut  also  on  as 
Indicated  In  Eqs.  12  and  27(b).  [Unfortunately,  an  explicit  relationship 
among  A^,  ply  elastic  constants,  lamination  parameters,  geometric,  and  loading 
variables  needs  to  be  determined  numerically  because  of  the  complex  algebra 
involved.]  We  remark,  that  the  maximum  Kj  and  Gj  appear  at  9  =  32.5°  and  the 
maximum  and  G2  occur  at  9  =  35°  in  the  [(±9)/(±9)/90°/90° ]s  graphite-epoxy 
system.  The  selection  of  [(±309)/(±30®)/90®/90° ]s  lamination  in  the  recently 
proposed  interlaminar  fracture  toughness  tests  [16,17]  using  edge-delaminatlon 
composite  specimens  is  therefore  a  proper  choice. 

5.2  Influence  of  Ply  Thlckners 

The  ply  thickness  is  an  important  geometric  variable  in  studying  the 
delamlnation  behavior  of  composites.  Changing  the  laminar  thickness  in  a 
laminate  alters  the  lateral  constraint  of  adjacent  plies  in  the  thickness 
direction  and,  thus,  directly  affects  the  interlaminar  crack  behavior.  For 
illustration,  fracture  mechanics  solutions  for  the  [450/-45°/-45°/45"] 
graphite-epoxy  composites  with  several  ply  thicknesses  are  presented  in  this 
section.  For  simplicity  and  without  introducing  further  complications,  the 
laminate  width  2b  and  thickness  2W  as  well  as  the  delamination  length  a  are 
kept  constant  as  before,  while  variation  of  ply  thicknesses  hj  and  h2  (with  h^ 
+  1>2  *  W  »  0.5  inch)  is  considered. 

As  shown  in  Figs.  8-10,  altering  the  ply  thickness  hj/W  has  appreciable 
effects  on  the  fundamental  behavior  of  the  delamination.  For  Instance,  when  a 
thick  outer  45°  ply  is  used,  say,  hj/W  >  0.8,  the  entire  delamination  surface 
is  in  contact  with  the  other  (Fig.  10),  leading  to  a  negative  Kj  with  Gj  - 
0.  Delamlnation  fracture  in  this  situation  is  governed  by  the  tecring-mode 
stress  intensity  factor  Kjjj  and  the  associated  energy  release  rate  G3  (Figs. 

8  and  9).  When  equal  ply  thickness,  i.e.,  hj  -  h2»  is  used,  the  in-plane 
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shearing-mode  stress  Intensity  Kjj  and  the  energy  release  rate  G2  are  about 
zero,  and  Km  has  the  maximum  stress  intensification.  However,  as  the  outer 
45s  ply  becomes  thin,  e.g.,  h}  <  0.2W,  the  surfaces  of  delamination  are  opened 
up,  i.e.,  c/a  *  0,  as  shown  in  Fig.  10.  In  this  situation,  contributions  of 
Kj  and  Kjj  to  the  local  failure  become  increasingly  appreciable;  Kjj  and  KtII 
are  of  the  same  order  of  magnitude  but  G2  and  Gj  are  still  relatively  small  as 
compared  to  G3.  In  the  limiting  case  of  h^/W  -*■  0,  the  in-plane  shear  could 
become  dominant  to  govern  the  delamination. 

5.3  Influence  of  Interlaminar  Crack  length 

For  the  symmetric  composite  laminates  with  either  [±0]s  or 
[<±9)/(±0)/9Oo/9O°]s  fiber  orientations,  the  interlaminar  crack-tip  behavior 
is  affected  by  the  size  of  delamination.  As  discussed  in  Section  5.1,  the 
delamination  in  the  [0/-0/-9/0J  graphite-epoxy  system  has  a  global  crack 
closure.  Figure  11  reveals  general  characteristics  of  the  crack  closure  in 
the  delaminated  [ 45° / — 45° /— 45®/ 45° ]  graphite-epoxy.  Under  a  uniform  axial 
strain  eQ,  the  crack  closure  increases  raonotonically  with  crack  length  until 
Che  delamination  becomes  about  two-ply  thicknesses.  As  the  delamination 
extends  further,  the  closure  length  approaches  an  asymptotic  value  of  c/h  * 
0.375,  indicating  that  crack  growth  is  governed  by  the  interlaminar  shear  ry2 
and  the  tearing-mode  stress  intensity  factor  Kjjj,  as  shown  in  Fig.  4.  The 
values  of  t^,  Kjj  and  G2  are  orders  of  magnitude  smaller  than  those  of  Tyz, 
Kjjj  and  G3  (Figs.  4  and  5).  Note  again  that  crack  closure  occurs  for  all 
a! h's  in  the  [45®/-45°/-45,/45®]  graphite-epoxy,  leading  to  an  identically 
zero  Gj  and  negative  opening-mode  stress  Intensity  factor  Kj  <  0.  In  the  case 
of  a  very  small  delamination  emanating  from  a  free  edge  under  rising  load,  Che 
interlaminar  crack  is  inherently  unstable  and  extends  rapidly  to  about  one  or 
two-ply  thicknesses  (i.e.,  the  plateaus  in  Figs.  4  and  5)  before  stable  growth 


2 


would  occur.  This  phenomenon  has  been  Indeed  observed  in  laboratory 
experiments  of  edge-delamlnatlon  tests  of  ((±30*]/±(30*)/90,79asl8  T300-5208 
graphite-epoxy  composite  laminates  (16). 

In  the  [(±30®)/(±30*)/90°/904]g  graphite-epoxy  containing  delamlnatlons 
between  -30°  and  90s  piles,  the  delaminatlon  Is  found  to  be  open  all  along  the 
crack  surface  for  any  given  a/h,  as  discussed  In  Section  5.1.  The  solutions 
for  Kj  and  as  a  function  of  delaminatlon  length  a/ b  are  given  In  Figs.  12 

and  13.  Also  shown  In  Fig.  13  are  the  results  reported  In  [16]  by  using  a 
conventional  nonsingular  finite  element  method  In  conjunction  with  a  virtual 
crack-closure  scheme  for  calculating  the  strain  energy  rates  [18].  The 
difference  Is  very  appreciable  to  warrant  the  necessity  of  using  the  advanced 
analytical  and  numerical  technique  for  solving  the  delaminatlon  problem 
accurately.  We  observe  that  local  delaminatlon  growth  In  the 
[(*30°)/(±30®)/90°/908]g  graphite-epoxy  Is  governed  by  an  Inherently  three- 
dimensional  mixed-mode  fracture  process  because  the  simultaneous  presence  of 
significantly  high  values  of  Kp  and  Kjh  shown  In  Fig.  12.  Owing  to  the 
vanishingly  small  A3,  the  value  of  G3  Is  found  to  be  negligible  In  comparison 
with  G}  and  G2  (Fig.  13).  We  remark  that  In  the  limiting  case  of  a  very  small 
delaminatlon  crack  (l.e.,  a/ h  *  0),  Interactions  occur  between  the 
delaminatlon  crack  tip  and  the  laminate  edge,  leading  to  a  slightly  higher 
value  of  Kjjj  as  shown  In  Fig.  12. 

A  salient  feature  shown  in  Figs.  4,  5  and  12,  13  is  that  the  stress 
Intensity  factors  Kj  and  energy  release  rates  Gt  become  Independent  of  the 
crack  length  <a/h  or  a/ b  as  the  delaminatlon  extends  beyond  a  few  ply 
thicknesses.  This  unique  feature  has  been  observed  and  used  in  the 
experimental  study  of  interlaminar  fracture  toughness  of  composite  laminates 
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[16,171.  Provided  thee  the  *  and  o,  are  acc.rat.ly  d.t.rdned  b,  appropriate 
aethode  and  a.  rbe  pra.ant  one,  the  air...  lnten.lt,  Unto,,  and  enetg, 
teleaae  rare.  ....dated  dtb  obaerved  edge  del.dn.tlon  Initiation  and  grovth 
be  naefnl  lor  ebaracterlalng  del.dn.tlon  Iradnt.  -nd  l.t.rlado.r  crack 

resistance  of  composite  laminates. 
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6.  SUMMARY  AMD  OONCLUSIOHS 

General  expressions  of  stress  and  deformation  fields  in  a  fiber  composite 
laminate  with  delaminatlons  have  been  obtained  from  presently  developed 
anisotropic  laminate  elasticity  theory  and  basic  concepts  of  Interlaminar 
fracture  mechanics.  Fracture  mechanics  parameters  such  as  stress  intensity 
factors  and  energy  release  rates  for  delaminations  with  different  crack-tip 
deformation  configurations  are  defined  in  a  manner  consistent  with  those  for  a 
homogeneous  crack  and  for  the  refined  model  of  an  interface  crack  between 
dissimilar  isotropic  solids  [7,8].  In  finite-dimensional  composite  laminates 
with  complex  lamination  and  geometric  variables,  an  advanced  numerical  method 
employing  special  singular  crack-tip  finite  elements  Is  developed  for  modeling 
the  delamination.  Exact  delamination  stress  singularities  obtained  from  the 
laminate  elasticity  solution  are  included  in  the  crack-tip  element 
formulation.  Solution  convergence  and  accuracy  have  been  studied  to  ensure 
the  validity  of  the  results  and  to  demonstrate  the  efficiency  and 
effectiveness  of  the  method.  To  illustrate  the  fundamental  nature  of 
composite  del imination,  numerical  results  are  shown  for  the  [0/— 6/— 0/0 ]  and 
[(±0)/(±0)/9O*/9O*]s  graphite-epoxy  laminate  systems  containing  edge 
delaminatlons  under  uniform  axial  extension.  Fracture  mechanics  parameters 
and  failure  modes  associated  with  the  composite  delamination  are  determined 
for  each  case.  The  basic  mechanics  and  mechanisms  of  delamination  are  studied 
for  the  composites  with  different  lamination  and  geometric  variables  and  crack 
parameters.  Based  on  the  results  obtained,  the  following  conclusions  may  be 


reached: 


(1)  To  study  the  fundamental  mechanics  of  delamination  la  a  composite 
laminate  having  a  large  number  of  plies  with  different  fiber  orientations, 
advanced  numerical  methods  such  as  the  present  one  are  essential  to  take  into 
account  the  crack-tip  singularity  and  the  large  number  of  lamination, 
material,  and  geometric  variables  involved. 

(2)  In  the  crack-tip  element  fonaulatlon,  the  Inclusion  of  delamination 
stress  singularities  determined  from  the  laminate  elasticity  solution  by  using 
an  eigenfunction  expansion  method  leads  to  very  accurate  solutions  with  a 
rapid  rate  of  convergence.  This  is  particularly  advantageous  for  studying 
delamination  problems  in  composite  laminates  with  complex  lamination, 
geometric,  and  material  variables. 

(3)  Since  the  singular  eigenvalues  and  their  algebraic  multiplicity  ar<* 
different  for  different  crack  tip  deformation  configurations,  the  crack 
element  formulation  and  solution  strategy  for  each  delamination  problem  need 
to  be  treated  on  an  individual  basis.  Stress  Intensity  solutions  and  energy 
release  rates  in  each  case  should  be  evaluated  in  accordance  with  the 
appropriate  delamination  models  given  la  Section  2. 

(4)  The  state  of  stress  and  deformation  in  the  vicinity  of  a  delamination 
crack  tip  are  three  dimensional  in  general.  The  asymptotic  solutions  can  not 
be  determined  accurately  by  using  classical  lamination  theory  nor  by  any 
approximate  methods  without  including  interlaminar  stresses  and  the  correct 
stress  singularities  associated  with  the  delamination.  The  current  laminate 
elasticity  solution  and  associated  numerical  method  provide  accurate 
information  on  the  singular  nature  of  the  crack  tip  and  complete  field 
solutions  for  the  delamination  problem. 
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(3)  The  slnguLar  de lamination  stress  field  nay  be  characterized  by  the 
presently  introduced  interlaminar  fracture  mechanics  parameters.  Craclt-tip 
stress  intensity  factors  and  energy  release  rates  G^  can  be  determined  only 
by  solving  the  complete  boundary-value  problem.  In  contrast  to  the  singular 
eigenvalues  which  depend  only  on  local  boundary  conditions  and  material 
properties  of  adjacent  plies,  the  and  are  functions  of  all  lamination 
and  geometric  parameters,  remote  loading  conditions,  and  crack  variables. 

(6)  Owing  to  the  local  nature  of  interlaminar  fracture,  the  mechanics  and 
mechanisms  of  delamination  growth  are  governed  by  the  crack-tip  stress 
intensity  factors.  In  angle-ply  [8/-0/-O/9]  graphite-epoxy,  crack-tip  closure 
occurs  and,  thus,  Kj  <  0;  delamination  growth  may  be,  therefore,  more 
intimately  related  to  interlaminar  shear  stresses  and  interface  shear  strength 
than  the  transverse  normal  stress.  In  the  cases  of  l(±8)/(±8)/90“/90*)g 
graphite-epoxy,  crack  surfaces  are  open  and  delamination  growth  is  controlled 
by  all  of  the  three-dimensional,  mixed-mode  stress  Intensity  factors,  Kj,  Kj-j, 
and  KIX1. 

(7)  While  the  govern  local  deformation  and  delamination  fracture,  the 
G^  are  related  to  the  global  structural  response  and  less  sensitive  to  the 
local  deformation  and  fracture.  For  example,  the  total  G  and  G-j  differ  only 
slightly  between  the  cases  of  a  delamination  with  a  closed  and  an  open  crack 
tip  in  the  1 0/— 9/— 0/9 J  graphite-epoxy  composite.  Thus,  Kj  may  provide  a  more 
sensitive  measure  and,  thus,  better  fracture  parameters  than  G1  for  evaluating 
the  composite  delamination  growth. 

(8)  Influences  of  lamination  variables  such  as  fiber  orientation  and  ply 
thickness  on  the  delamination  behavior  are  significant.  Changing  fiber 
orientation  generally  alters  failure  modes  appreciably.  For  instance,  in 
l(*fl)/(ifl)/90*/90#l9  graphite-epoxy,  the  opening-mode  (Kj)  dominated 


delaainatlon  fracture  changes  into  a  shearing-node  (Kjj)  controlled  failure, 
as  the  9  be  coses  greater  than  40*.  Also,  increasing  ply  thickness  hj/W  in 
[9/-0/-0/9]  graphite-epoxy  affects  the  failure  node  from  a  shearing-dominated, 
open  delamination  fracture  to  a  closed  one. 

(9)  Stress  Intensity  solutions  and  energy  release  rates  appear  to  be 
independent  of  crack  length  during  the  growth  of  an  edge  delaelnation ,  as  long 
as  the  crack  exceeds  a  length  of  few  ply  thicknesses  [3,19].  This  unique 
feature  in  composite  edge-delaalnatlon  is  being  used  for  evaluating  of 
lnterlaalnar  fracture  toughness  of  fiber  composite  laminates  [16,17],  provided 
that  the  and  can  be  calculated  accurately  by  using  advanced  analytical 
methods  such  as  the  present  one. 
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TABLE  1 

Delamination  Stress  Intensity  Factors  K^+  Associated  with 
Finite-Element  Mesh  Refinement  in  I45*/-45*/-45*/45*]  Graphite-Epoxy 


No.  of 
Elements 

Element  Size 
(Ah/h) 

* 

KI 

KII 

Km 

72 

0.25 

-0.6458E-1 

0.3424E-2 

-0.7114E  0 

84 

0.125 

-0.6219E-1 

0.2463E-2 

-0.6967E  0 

96 

0.0625 

-0.6141E-1 

0.1796E-2 

-0.6894E  0 

108 

0.03125 

-0.6152E-1 

0. 1401E-2 

-0.6855E  0 

+Kw  are  normalized  by  106e  (psi-/in)  . 

o 

* 

Closed  crack  tip  with  c/a  ■  0.34. 


e 
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TABLE  2 

De lamination  Energy  Release  Rates  C^f  Associated  with 
Fialte-Eleaeoc  Hash  Ref ioeaenC  in  (45*/-45*/-45*/45* )  Graphite-Epoxy 


No.  of 

Elements 

Element  Size 
( Ah/h) 

ft 

C1 

C2 

C3 

72 

0.25 

0. 

0.6855E-6 

0.1533E  0 

84 

0.125 

0. 

0.5528E-6 

0.1519E  0 

96 

0.0625 

0. 

0.4619E-6 

0.1512E  0 

108 

0.03125 

0. 

0.3943E-6 

0.1508E  0 

are  normalized  by  10*c2  (lb- in/in2). 
Closed  crack  dp  with  c/a  *0.34 
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TABLE  3 

* 

Influence  of  Fiber  Orientation  on  Crack-Tip  Cloture  Lengthf 
for  a  Da lamina t ion  in  f ©/-«/-©/ 0 }  Craphlte-Epoxy  Composite 


9 

c/a 

0* 

0.0 

15* 

0.2202 

30* 

0.2834 

45* 

0.3402 

60* 

0.3603 

*h  •  h_  •  0.25  in.,  b  »  2  in.; 
a  ■  0.25  in. 

*«t  strain  c  •  10~6 
o 
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TABLE  4 


Effects  of  Fiber  Orientations  on  Stress  Intensity  Solutions 
for  Delamination  in  [0/-9/-Q/9]  Graphite-Epoxy  Composite* 


0 

Ki 

KII 

KIII 

0° 

0. 

0. 

0. 

15° 

-0.3782E-1 

0.3866E-2 

-0.2379E  1 

30° 

-0.1072E  0 

0.8744E-2 

-0.1943E  1 

45° 

-0.6152E  0 

0 . 1401E-2 

-0.6855E  0 

ON 

o 

0 

-0 . 8599E-2 

0.2034E-3 

-0.9270E-1 

+hl  =  h2  =  °-25  in.,  b  =  2.0  in.;  a  =  0.25  in. 
are  scaled  by  106eQ  (psi-ZErTl. 
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TABLE  5 

★ 

Effects  of  Fiber  Orientation  on  Strain  Energy  Release  Rates 
for  a  Delamination  in  [8/-8/-0/0]  Graphite-Epoxy  Composite' 


8 

G2 

G3 

0° 

0. 

0.0 

0.0 

15° 

0. 

0.6059E-5 

0.2111E  1 

0 

O 

fO 

0. 

0 . 2748E-4 

0.1074E  1 

45° 

0. 

0.3943E-6 

0.1508E  0 

60° 

0. 

0.1425E-7 

0.3597E-2 

\ 

\ 

=  h2  =  0.25  in.,  b  =  2  in.;  a  •  0.25 
are  scaled  by  (106e^) [lb-in/in2] . 

in. 

is  identically 

zero  due  to  crack-tip 

closure. 
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m  hn  ■  0.0054  in.,  b  -  0.75  in.). 
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(b) 


Fig.  I  (a)  Six-Node  Quasi  Three-Dimensional  Crack-Tip  Element  with 
Singular  Derivatives;  (b)  Eight-Node  Quasi  Three-Dimensional 
(nonsingular)  Isoparametric  Element. 


Arrangement  of  Special  Crack-Tip  Elements  for  the  Composite 

n _ L 1 _  1 


Fig.  3  Overall  and  local  Crack-Tip  Fini te-Element  Meeh  Arrangement*  for 
Modeling  Oelami nation  in  {43V-43V-4 5*/43*J  Graphite-Epoxy 
Composite. 


Singular  F.E.M. 

(Partially  Closed  Crack) 
Eigenfunction  Expansion 
(Fully  Opened  Crack) 


Fig.  4  Comparison  of  Stress  Intensity  Solutions  Obtained  by  Two  Different  Models  for  Delamination 
in  (A50/-A5°/-A5,,/A5o ]  Graphite-Epoxy  Subjected  to  Uniform  Axial  Strain  e  ■  eQ  (hi  ■  hi  ■ 
h  -  0.25  in.,  b/h  -  8  in. ,  a  -  0.25  in.) . 


Singular  F.E.M.  (Partially  Closed  Crack) 
Eigenfunction  Expansion 


Comparison  of  Energy  Release  Rates  Determined  by  Two  Different  Models  for  Delamlnatlon  In 
|45°/-450/-A50/45° ]  Graphite-Epoxy  Subjected  to  Uniform  Axial  Strain  e  ■  e  (hi  *  ho  *  h 


K,  /{IO*€0)  (psi-Vln) 


Fig.  8  Influence  of  Ply  Thickness  hi/W  on  Scree*  Intensity  Factors  for 
Delaalnation  in  (45*/-45*/-45*/45* J  Graphite-Epoxy  Subjected  to 
Onifotfa  Axisl  Strain  ez  •  t0  (hj  +  h2  -  W  •  0.5  in.,  2b/W  -  8,  a  - 
0.25  in.). 


X 


CiUfl  * 

OF  POOR  QUALITY 


h,/W 


Pig.  9  Influence  of  Ply  Thickness  h>/V  on  Energy  Release  Rates  Ct  for 
DtUalnatlon  In  (*5*/-45*/-45*/*5* J  Graphite-Epoxy  Subjected  to 
Onlfora  Aatlal  Serein,  c.  -  ta  (hi  ♦  h2  -  V  -  0.5  In.,  2b/V  -  8,  3  - 
0.25  In.). 


c/a 


[45°/-  45%  45°/45°]  Grophite  -  Epoxy 


FI*.  10  Crack-Tip  Cloture  Length  as  a  Function  of  Ply  Thickness  hj/V  for 
Oelaal  nation  In  (45V-*5*/-*5*/45*)  Graphite-Epoxy  Subjected  to 
Onlfora  Axial  Strain  e,  -  tQ  (hj  ♦  h2  -  W  -  O.S  In.,  2b/V  -  8,  a - 
0.25  In.). 


Crack-Tip  Closure  Length  c/h  as  a  Function  of  Delamlnat ion  Size  a! h 
in  (45°/-450/-A50/450 1  Graphite-Epoxy  Subjected  to  Uniform  Axial 
Strain  e_  ■  e  (hi  »  ho  =»  h  »  0.25  in.,  b/h  -  8,  a  -  0.25  in.). 


[±30o/±30°/90o/90°]s  Graphite -Epoxy 
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(ui^-!Sd)  (°3,0I)/!M 


I 


ess  Intensity  Solutions  with  Delamination  Length 
)/(±30°)/90o/90o"]s  Graphite-Epoxy  Subjected  to  Uniform 
c_  -  e„  (hi  -  ho  ■  ...  ■  hii  “  0.0054  In.,  h  "  0,75 


Present  Singular  F.E.M. 
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OF  POOR.  QUALITY 


(.uiyui-qi)  (?3,oi)/!9 


* 


Change  of  Strain  Energy  Release  Rates  G^  with  Delamination  Length 
alb  in  [(±30<’)/(±30°)/90®/90ols  Graphite-Epoxy  Subjected  to  Uniform 
Axial  Strain  e_  -  tn  (hi  -  ho  “  ...  ■  h»i  -  0.0054  in.,  b  »  0.75 
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APFESDIX  1 

Shape  functions  Mg(p,C;6)  for  Che  Singular  Del ami nation  Crack-Tip  Element 

Corresponding  to  a  singular  eigenvalue  6,  the  stiape  functions  N  in  local 

«  s 

triangular  polar  coordinates  (p.O  for  the  six-node  conforming  crack-tip 
element  in  the  delamination  problem  are  shown  [4]  to  have  the  following 


expressions: 

Nu  -  1  +  ((2-2_,5)p-p5+1l/(2"6-l),  (Al-1) 

Nl(3+i)  *  2"4p(1-C)/(2_<5-1)  -  p5+1[2e(l-0+(l-0/(2"6-l)],  (Al-2) 

Nl(6+i)  ’  2“4pC/(2"6-1)  -  p5+1{2C(l-0+C/(2"6-i)I,  (Al-3) 

Ni(9+i)  *  "  [2p(1-5)+2p6+1(1-5))/(2“<S-1),  (Al-4  . 

Ni(l2+i)  ”  4P6+1?(1-«)*  (Al-5) 

N1(i5+i)  -  -  2[p£-p6+1£l/(2~6-l>,  (i  -  1,2,3),  (Al-6) 


where  the  singular  eigenvalues  6’s  are  determined  in  accordance  with  the  local 
crack-surface  boundary  conditions  discussed  in  Section  2.1,  and  the  p  and  £ 
are  related  to  the  global  coordinates  by 

(x2~x  j)  tan  4  -(y2_yp 

5  ”  (y3~y2^  tan 

P  -  r/f(£). 


(Al-7) 

(Al-8) 


in  which 


f(£)  -  ((x2-x1)2+(y2-yl)2+2U(x2-xl)(x3-x2)+(y2-y1)(y3-y2)] 


+£2l(x3-x2)2+(y3-y2)2]}1/2. 


(Al-9) 
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APPSMD1X  2 

Iteration  Scheme  for  Partially  Cloned  De  land,  nation 

For  a  delaainatlon  with  crack  surfaces  In  finite-length  contact,  the 
following  iteration  scheme  is  used  for  determining  contact  stress  and  crack 
closure  length: 

(i)  Assume  an  initial  contact  length  cj,  and  solve  Eq.  21  with  the  side 
conditions  Eqs.  23(a)  and  23(b)  for  contact  stress  and  displacement. 

(li)  Check  the  solution  against  preset  criterion  nf.  If  the  compressive 
stress  F^  -  <J2(-ci,*)  gives  a  nj  In*  =  I (Fi“F0)/Fo I >  and  Fo  ls  3 
properly  selected  sealing  factor,  for  example,  FQ  ls  set  as  10^  eQ 
(psi),  in  the  present  calculation.]  such  that  ni  <  Hf  and  if  the 
displacement  field  ls  admissible  (i.e.,  no  overlapping  or 
interpenetration  beyond  the  contact  region),  we  set  the  crack  closure 
length  c  “  Cj  and  terminate  the  iteration. 

(ill)  If  the  two  constraint  conditions  are  not  met,  a  new  contact  length  C2  > 

cj  is  assumed  and  the  procedure  of  (i)  and  (li)  is  repeated  for  F2  and 

u2. 

(lv)  If  F^  <  0  (i  *»  1,2)  and  |n2|  >  |  HjJ  >  Hj,  the  next  assumed  length  C3 
should  be  C3  <  cj,  and  repeat  (iii)  and  (iv). 

(v)  If  either  F^  <  0  and  |nj|  >  |n2|  >  \  or  FjF2  <  0  and  |nt|  >  e^,  the 

next  assumed  contact  length  ls  set  as 

c3  -  c2  +  F2(c2-c1)/(Fl-F2).  (A2-1) 

The  iteration  from  (i)-(lil?  continues  until  sufficient  accuracy  ls 


reached 
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(vl)  As  Che  constraint  conditions  are  satisfied  and  the  difference  between 
assuaed  contact  lengths  Is  Ic^-Cnj  <  £c,  the  Iteration  is  terainated 
withe  -  V2  lcn+l+cj* 


